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ABSTRACT 

We have developed a new method for post-Newtonian, high-precision integration of 
stellar systems containing a super-massive black hole (SMBH), splitting the forces on 
a particle between a dominant central force and perturbations. We used this method 
to perform fully collisional iV-body simulations of inspiralling intermediate-mass black 
holes (IMBHs) in the centre of the Milky Way. We considered stellar cusps of different 
power-law indices and analysed the effects of IMBHs of different masses, all starting 
from circular orbits at an initial distance of 0.1 pc. 

Our simulations show how IMBHs deplete the central cusp of stars, leaving be- 
hind a flatter cusp with slope consistent with what has recently been observed. If an 
additional IMBH spirals into such a flat cusp, it can take 50Myr or longer to merge 
with the central SMBH, thus allowing for direct observation in the near future. The 
final merger of the two black holes involves gravitational wave radiation which may 
be observable with planned gravitational wave detectors. 

Furthermore, our simulations reveal detailed properties of the hypervelocity stars 
(HVSs) created, and how generations of HVSs can be used to trace IMBHs in the 
Galactic Centre. We find that significant rotation of HVSs (which would be evidence 
for an IMBH) can only be expected among very fast stars (v > 1000 km s -1 ). Also, 
the probability of creating a hypervelocity binary star is found to be very small. 
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1 INTRODUCTION 

In recent years, observations of the Galactic Centre revealed 
a number of interesting phenomena. Analysis of the motions 
of stars within 1.2 arcsec of the Galactic Centre, for exam- 
ple, showed that the Milky Way contains a super-massive 
black hole (SMBH) of mass « 3.5 x 10 6 M B , surrounded 
by an extremely dense stellar cluster dSchodel et al.l [20031: 
iGhez et alj|2005h . 

Since mo st massive galaxies are now believed to hos t 
such SMBHs jMagorrian et al.ll 19981 : [Cebhardt et alj|2000h . 
due to its proximity (~ 8kpc), the centre of the Milky Way 
provides a unique environment to study the dynamics of 
SMBHs and the surrounding stellar systems. 

While SMBHs have been directly observed as radio 
sources and compact massive objects, the existence of black 
holes of masses 1O 2 M < M BH < 10 4 Af Q , so-called 
intermediate-mass black holes (IMBHs), is still a matter of 



debate (e.g. Hopman. Portegies Zwart fc Alexander! [20041 : 
iMiller fc ColbertjhoOll ; IPortegies Zwart et al.ll2004h . 

IMBHs could form in galaxies through runaway col- 
lisions of stars in star clusters, giving rise to ultra- 
luminous X-ray sources llGiirkan. Freitag fc Rasiol 120041 : 
IPortegies Zwart et al.ll2004l : iBaumgardt et al.ll2006bl '). They 
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could later be brought into the centres of galax- 

ies through dynamic al friction of the star clusters 

|Portegies Zwart et alj |2006| ). Some of the dynamical pro- 
cesses at the Galactic Centre may be well explained by the 
presence of an IMBH, such as the existence and kinematics 
of yo ung stars in the central 0.1 pc l|Hansen fc Milosavlievid 
l2005t) . Mergers of IMBHs with the central SMBH would 
be an important source of gravitational waves, detectable 
with the next generation of gravitational wave detectors, 
like e.g. the Laser Interferometer Space Antenna (LISA, 
iBender et ai1ll99Sl ). 

During the inspiral process of an IMBH, close encoun- 
ters with stars from the central cluster eject some of them 
with velocities large enough to leave the gravitational poten- 
tial of the Milky Way, so-called hypervelocity stars (HVSs) . 
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iHillsl (|l988l ) was the first to show that the ejection of HVSs is 
a natural consequence of galaxies hosting SMBHs. Recently, 
several HVSs have been discovered in the Galactic ha lo 
llBrown et al.ll2005l ; iHirsch et ai]|2005l ; iBrown et aDl2006al lbl. 
l2007al lbl). Except for one star which might have been 



eject ed from the Large Magellanic Cloud ( Edel mann et al.l 
2005), the travel times of all HVSs are short enough 
that the stars could have been ejected from the Galac- 
tic Ce ntre within their lifetimes, confirming Hills ' predic- 



tions. 
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and iLevin 
ejected almost 
railing IMBHs, explaining the near-isotrop i c distr ibution 
of th e so far detected HVSs dBrown et all l2007al ). How- 
ever, ISesana. Haardt fc Madaul ^Ool )" found that the rel- 
atively small medium velocity of the (admittedly limited) 
sample of HVSs detected so far favours other eje ction 
mechanisms like the t i dal breakup of stellar binaries ( Hills 
19881; lYu fc Tremaine, 120031 ; IPerets. Hopman fc Alexander 



2007) or the scattering of stars by stellar-mass black holes 



i Miralda-Escude fc Gouldll2000l ; IO'Learv fc Loebll2008l ). 

In this work, we present new results on the inspiral of 
IMBHs, using a novel method for studying the dynamics 
of Galactic Centre-like systems where stars orbit a central 
SMBH on weakly perturbed Keplerian orbits. In Section [2] 
we describe the details of our integration method. In Sec- 
tion [3] we describe our set of N-body runs with IMBHs in 
the Galactic Centre. Our results are given in Section [3] and 
discussed in Section [5] 



2 INTEGRATOR FOR SMBH-DOMINATED 
SYSTEMS 

2.1 Requirements 

Stellar systems around SMBHs are a unique environment: 
Dynamically, they behave like the solar system in that stars 
move along weakly perturbed Keplerian orbits around a 
massive central body. Integration of planetary systems is 
usually done by symplectic integrators, which yield a very 
good energy conservation for nearly circular orbits but fail 
for eccentric orbits with a wide range of semimajor axes. 

The structure of the stellar system around a massive 
black hole is also similar to a star cluster, since it comprises 
a large number of stars of similar mass with wide ranges of 
eccentricities and central distances. Good progress has been 
made in the development of high precision integrators for 
the dynamical simulati on of star clusters, usuall y based on 
the Hermite algorithm dMakino fc Aarsethll 19921 ). 

A drawback underlying Hermite integrators is their in- 
ability to differentiate massive objects dominating the sys- 
tem's dynamics from relatively small masses causing only 
small perturbations to the particles' orbits. 

Although integrators have been developed in both the 
areas of planetary systems and stellar clusters, a new ap- 
proach is therefore required for the Galactic Centre. 

We have developed an integrator bhint specialized for 
dynamical processes in the vicinity of a SMBH, making use 
of the information that the black hole dominates the motion 
of the stars. For this new mechanism, we retain the Her- 
mite scheme as a basis, providing very fast calculations on 



the GRAPE special-purpose hardware ( M akino et al.l [20031 ) 
which is essential for the large values of we use. 

2.2 Basic integration algorithm 

The basic idea for our new method is to split the force calcu- 
lation between the dominating central force (exerted by the 
SMBH) and the perturbing forces (due to the cluster stars). 
It is comparable to the so-called mixed v ariable symplectic 
methods (MVS, I Wisdom fc Holmanlll99lh . as it makes use 
of Kepler's equation to integrate along the orbit. However, 
our method is not symplectic, as it is based on the Hermite 
scheme to allow for use of the GRAPE special-purpose hard- 
ware. Furthermore, we assume the dominating SMBH to rest 
at the centre and hence do not need Jacobi coordinates but 
calculate the orbital motion in Cartesian coordinates. 

A different a dvancement of th e idea of MVS was also 
made recently bv lFuiii et al.l (|2007l ) for simulation of a star 
cluster and its parent galaxy, separating forces between star 
cluster particles and galaxy particles. 

The underlying Hamiltonian for a Newtonian problem 
of N particles with masses mo..mn_i, positions ro..rM-i 
and velocities vq-.vm-i is 



H — \P m i v i _|_ (jTn i m j 



(1) 



where rt 



-Vi, 



\rij\, and G is the gravitational 



constant. As the central object is much more massive than 
the orbiting particles, we will assume it does not move with 
respect to the system's centre of mass (r () = 0)Q In this 
case, the total force acting on a particle can be split into 
the unperturbed motion around the central object, and the 
perturbing force, as 



-Gm ^-J2 Gva^ 



(2) 



for i = 1...N- 1. 

The new position of a particle after a time-step At can 
be expressed as 



1 ,. A 2 

r + rAt +—r K Ar 



1 „(3) 

3! 

„(3) 



At 6 + . 



(3) 



The first line of equation ([3]) is solved directly, assum- 
ing an unperturbed motion along a Keplerian orbit and 
making use of Kepler's equation. The second line is calcu- 
lated using the Hermite i ntegration scheme as described in 
iMakino fc Aarsethl (|l992l ). To adapt this scheme to our in- 
tegrator, the predictor step for the particles to be moved 



For our models, the expected motion of the SMBH relative to 
the system's centre of mass is of the order of 3 X 10 — 4 pc. Only 
few of the stars have a pericentre distance within that range, and 
most of these will move around with the SMBH, as their orbital 
period is much shorter than the expected S MBH motion time- 
scale (cf. Chattcriec, Hernquist & Loeb 2002). Hence, the motion 
of the SMBH can be neglected. 
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needs to be changed to 
»Ved = ri,K + \a P At 2 



0.001 



-apAt 3 
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(4) 
(5) 



1 a 1 . A 2 

« P red = »i,k + -apAt+ -apAt 

where ri.jc and Vi,k are the results of Kepler integration. 

To determine the size of a particle's next time-step, the 
following criteria are used: 

(i) The particle's central motion. The accuracy of the Ke- 
pler orbit integration does not depend on the step size, how- 
ever the faster change of the particle's velocity at perihelion 
requires a smaller time-step as it causes faster changes in 
perturbing forces. Hence we introduce a scale-independent 
upper time-step limit as 



Atkep = 



2tt 



(6) 



iVstep V GA/sMBH 

where r is the particle's central distance, Msmbh is the 
SMBH mass and JV at e P a constant which we set equal to 
50. For a circular orbit, N B t ep gives a lower bound to the 
number of steps per orbit (more steps are needed at the 
pericentre of eccentric orbits). 

(ii) The relative change of pe rturbing forces . Here, we ap- 
ply the criterion introduced by lAarsethl (|l985l ) , 



. i |a||a| + |a|- ,_. 

Atpcrt = \j n \a\\a^\ + \ar () 

and use 77 = 0.5. lAarsethl |l985h suggests a value of r\ = 0.02; 
however, in our case changes of perturbations are usually of 
rather little effect on a particle's trajectory, as the orbital 
motion around the SMBH is dominating. Thus to retain a 
certain level of energy conservation, the time-steps can be 
much larger compared to a standard Hermite scheme, as the 
dominating orbital energy is conserved in any case by our 
method. 

(iii) To prevent jumps in energy by not handling close en- 
counters correctly, we need to consider changes in the par- 
ticle's vicinity (i.e. approaches and close encounters). We 
assign a neighbour sphere to each particle and check for fast 
approaches and close encounters before we perform a step. 
Thus, we prevent sudden changes of perturbing force deriva- 
tives. A particle's close encounter step size At cnc is set to the 
maximum of all step sizes, during which the particle will not 
approach any other particle by more than a factor of ~ 3. 

The definitive time-step of a particle i is defined as the 
minimum of Aikcp, At pclt , and At cnc , reduced to the next 
smaller integer power of two 



Ati = 2 



«)J 



(8) 



In our simulations, all particles are treated as point 
masses. We do not consider tidal disruption during our cal- 
culations, but track close encounters between two particles. 

Our simulations show that binary formation in a 
SMBH-dominated cusp can be neglected, therefore we do 
not use any regularization methods. 

2.3 Comparison with other integration methods 

High-precision symplectic integrators are available for differ- 
ent kinds of problems, such as the time-symmetric adaptive 
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Figure 1. Comparison of our new integrator method BHINT (solid 
lines) with the standard Hermite scheme (dotted lines) for a clus- 
ter of 1000 stars with masses 3 — 10 Mq on orbits with a range of 
eccentricities and semimajor axes between 0.05 and 0.1 pc around 
a SMBH of mass 3.5 X 10 6 Mq. For a given average number of 
steps per orbit, BHINT is almost a factor of 100 more accurate 
than the standard Hermite method. 



time-s tep mechanism introduced by Mikkola & Tanikawa 
1 19991) , which was advanced by iLevin. Wu fc Thommesl 
(2005) for use in a three-body problem to study the effects 
of an inspiralling IMBH in the Galactic Centre on a single 
star, or the indi vidual time-step multiple variable symplec- 
tic integrator bv lSaha fc Tremaind (| 19941 ). which is designed 
for planetary systems. However, symplectic integrators fail 
when dealing with Galactic Centre-like systems, as these 
comprise a wide range of semimajor axes and eccentricities, 
and no symplectic integrators with individual and adaptive 
time-steps have been found as yet. 

Therefore, we can only compare our new method to an 
integration scheme fo r large stellar clusters, a mong which 
the Hermite method l|Makino fc Aarsethlll992l) is the most 
widely used method for star cluster integration. 

Fig. [1] compares the performance of our new method 
bhint with a standard Hermite scheme. It shows the rela- 
tive energy error for integration of a cluster of 1000 stars 
with masses 3—10 Mq on orbits with a range of eccentric- 
ities and semimajor axes between 0.05 and 0.1 pc around a 
SMBH of mass 3.5 x 10 6 Mq. For a given average number of 
(100, 200,400) steps per orbit, one can see that our method 
is almost a factor of 100 more accurate than the standard 
Hermite method, as it does not accumulate an error in the 
orbital motion around the SMBH. 

The secular ene rgy error of the Hermite scheme ha s 
been investigated bv lKokubo. Yoshinaga fc~M akino ( 1998h . 
who also introduced a multiple-evaluation Hermite scheme, 
P(EC)™, to overcome this drawback for low-eccentricity sys- 
tems with rare close encounters (such as planetary systems). 
Our new method obviates the energy error of the orbital mo- 
tion (which is calculated exactly by solving Kepler's equa- 
tion), but only retains the energy error due to interactions 
between the orbiting particles, which in systems dominated 
by a central body is much lower than the total energy. 
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2.4 Post-Newtonian extension 

As an IMBH spirals into the Galactic Centre due to dy- 
namical friction, it can eventually reach a highly eccentric 
orbit with a relatively small semimajor axis, a regime where 
relativistic effects beco me important. For example, accord- 
ing to equation (5.10) of lPeterd (|l964T ). a 3 x 10 3 M Q IMBH 
orbiting a 3.5 x 10 6 M SMBH on a 10 3 pc orbit with eccen- 
tricity 0.97 will merge with the SMBH within ~ 1 Myr due 
to gravitational radiation. The actual value is even smaller, 
as the relativistic pericentre shift prolongates the pericentre 
passage. 

To account for relativistic effects, we extend our inte- 
gration method by post-Newtonian (VN) correction terms. 
This changes the calculated acceleration to 

a— do + c~ 2 a,2 + c~ 4 a4 + c~ a$ +0(c~ 6 ). (9) 

Newtonian lVAf 27? Af 2.5T>JV 

pericentre shift GW 

The actual formulae for the VN correction terms are taken 
from iBlanc hct (2006, equation (168)), simplified to our case 
of a SMBH at rest. To use these in our Hermite scheme 
for p erturbing forces, we also calculated the force deriva- 
tives. iKvrrjij^rriaro^Se^arie^^Srjmze^ |2006Tl were the first 
to include these post-Newtonian corrections into a standard 
Hermite integration method. 

While the 2.5VAf term corresponds to the emission of 
gravitational waves and thus does not conserve energy, the 
1VN and 2VN terms are responsible only for pericentre 
shift and do not change the semimajor axis and eccentric- 
ity. However, even these two terms lead to a large periodic 
change in the Newtonian energy. 

For these reaso ns, we track the integrator's energy con- 
servation similar to lAarsethl (|2007T ) by evaluating the path 
integral 

AE = mva Gw dt. (10) 



Run 


-WlMBH 


M* 


N 


a 


A 


10 3 


30 


19,470 


1.75 


B 


3 x 10 3 


30 


15,779 




C 


10 3 


30 


13,549 




D 


3 x 10 3 


2.5 - 8.5 


58,483 


1.4 



Table 1. Parameters used in N-body runs. Masses are given in 
M Q , N is the number of particles. The SMBH mass is 3 X 1O 6 M 
in all runs. Runs B and C take the results of runs A and B, 
respectively, as input, plus a follow-up IMBH. They thus start 
with an a = 1.75 density profile at large distances and a slightly 
(run B) or strongly (run C) depleted cusp at smaller distances. 

3 DESCRIPTION OF THE RUNS 

All runs were performed with the code descr ibed above on 
the GRAPE6 computers ijMakino et al.ll200St ) of Bonn Uni- 
versity. Our runs contained three different components: A 
central SMBH, an IMBH, and a cluster of stars [varying 
in size and stellar mass, but fo llowing the stellar mass es- 
timated by iGenzel et all (|2003h within the inner 0.4 pc]. In 
all simulations, the SMBH was at rest at the origin and 
had a mass of Msmbh = 3 x 1O 6 M0, similar to the mass 
of the SMBH at the Galactic Centre (|Schodel et all 120031 ; 
iGhez et all 120051 ). The mass of the IMBH was varied be- 
tween (1-3) x 10 3 Mq in the different runs. All IMBHs 
moved initially in circular orbits at a distance of 0.1 pc from 
the SMBH. Details are shown in Table [1] 

The stellar cluster was set up such as to be initially 
non-rotating and fulfill an isotropic velocity distribution, 

e*v 2 =e^ =E"HE" 2 d 2 ) 

(e.g., iBinnev fc Tremaind 1 19871 . equation (4-53b)), which 
translates into a uniform distribution of e 2 and a mean cen- 
tral distance r = 5/4 a (where e and a are a star's eccen- 
tricity and semimajor axis, respectively). The initial density 
profile of the stellar cluster was a power law with different 
exponents for the different runs (see Table [TJ. 



To account for relativistic effects, we replace the Keple- 
rian force r^x in equation ((2)1 by the relativistic Keplerian 
force 



-Gmo^-j+c 2 a,2 + c 4 a4 + c J a 5 +C(c 6 ). (11) 



Since no analytical solution for the post-Newtonian two- 
body problem is known, we integrate this equation us- 
ing a standard Hermite scheme. In contrast to the Kepler 
equation, this is not exact, and especially since the post- 
Newtonian forces change rather quickly, we replace every 
Kepler step by 256 steps of integrating equation pip . Only 
after performing these 256 steps, the perturbing forces are 
applied to a particle. 

As the integration of the post-Newtonian terms is time- 
consuming and we are mainly interested in the IMBH inspi- 
ral, we will only apply it to the IMBH, and only when its 
estimated merging time falls below 100 Myr. It is also possi- 
ble to apply this integration to all particles below a certain 
threshold, trading computational speed for higher accuracy 
of stars on close eccentric orbits. 



4 RESULTS 

4.1 Predicted IMBH inspiral 

Due to dynamical friction, an IMBH sinks into the centre 
of a stellar cusp. The frictional drag on the IMBHs can be 
estimated by [see IBinnev fc Tremaind l)l987l '). equation (7- 
18)]: 



dv 
~dt 



4tt In AG 2 p(r)M n 



„, v , 2X „ X 2 
ert(X) -=e 

\/7T 



(13) 



where p(r) is the background density of stars, In A the 
Coulomb logarithm and X = v/(\/2a) is the ratio between 
the velocity of the IMBH and the (ID) stellar velocity dis- 
persion <t, which can be inferred from the model data. For 
a stellar density profile oc r~ a with 1.4 < a < 1.75, one gets 
X « 1.2. 

For a power-law density profile p(r) — por~ a , assuming 
the IMBH moves on a circular orbit, and introducing the 
constant 



: 87rlnA 



GpoMlMBH 



erf (X) 



2X 
V 71 " 
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equation (|13[) can be rewritten as 



dv 

r = AiiMBH 77" 
dt 



= -yMlMBH^r 



3/2- 



(14) 



Since the rate of angular momentum change is equal to 
dL/dt = Fr /Mimbh and since the angular momentum itself 
is given by L — rv c = -v/GA/smbhT, we have 



Fr 



dL 



v c dr 



- = -VGM SMm r = --, 



Mimbh 

and inserting equation (|17|l leads to 
dr 



a-5/2 



dt 



A. 



(15) 



(16) 



Solving this equation with initial condition r(0) = ro 
gives the radius reached at time t: 



r(t) 



. Q -3/ 2 



.4 



(a -3/2)*)' 



(17) 



From this it can be seen that for cusps steeper than a — 
1.5, the inspiral process terminates in finite time, whereas 
in shallower cusps the IMBH will never reach r = 0. 



4.2 IMBH inspiral in a Bahcall-Wolf cusp around 
a SMBH 

Theoretical arguments and N-body simulations have shown 
that a stellar system around a SMBH evolves into a 

cusp with an a = 1.7 5 power-law density distribution 

(iBahcall fc WolJ Il97rj: IBaumgardt, Makino fc Ebisuzakil 
12004 al ibi: IPreto. Merrill ,v Spurzemll2004r i. 

Following these theoretical predictions, we started 
our calculations with a 10 3 Mq IMBH on a circular or- 
bit around a 3 x 10 6 Mq SMBH with semimajor axis 
ffliMBH = 0.1 pc (run A). Fig. [2] shows that the inspiral 
process follows the theoretical description (equation (|17[H 
very well. We find a best fit for all our runs if we as- 
sume a Coulomb loga rithm of In A = 7.5, which is in 
good agreement with iMatsubavashi. Makino fc Ebisuzakil 
(2007). The actual inspiral deviates from the prediction 
once the IMBH has started to deplete the inner cusp 
and the density profile no longer follows an r -1,75 law. 
These results agree very well with tho s e obta ined by 
IBaumgardt. Gualandris fc Portegies Zwartl (|2006al ). 

From this point, the IMBH's central distance de- 
creases only slowly, as two-body encounters are rare. 
Gravitational wave (GW) emission at the pericen- 
tre becomes important after 8 Myr, when the IMBH 
acquires a highly eccentr ic orbit due to i n terac- 
tions with passing stars JHc^ man fc Alexander] |2006; 



IMatsubavashi. Makino fc Ebisuzakil " 2007 1. This leads to 
further decline of the semimajor axis, and later on to a 
circularization of the orbit. The inspiral time due to GW 
radia tion is about 1 Myr (using equation (5.14) from lPeters! 
1 1964? ). but orbit circularization delays this process. 

After the IMBH in run A has started to deplete the 
inner cusp and reached a central distance of 5 x 10 -4 pc, 
we start a subsequent run (B), assuming that a follow-up 
3 x 10 3 Mq IMBH starts from a distance of 0.1 pc. Due to 
the higher mass, the inspiral process is much faster (Fig. 
f2|. After 1 Myr, the IMBH reaches the innermost region, 
where the density is decreased and the inspiral process is 
slowed. After 1.5 Myr, the IMBH acquires a highly eccentric 
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Figure 2. Evolution of the inspiralling IMBHs' eccentricities 
e lMBH an d semimajor axes aimbh ■ The solid and dashed lines de- 
pict the results of runs A (Mjmbh = 10 3 Mq) and B (AfiMBH = 
3 X 10 3 Mq), respectively, while the theoretical predictions are 
plotted in grey. Once the eccentricity increases, GW emission be- 
comes important at the pericentre, leading to circularization of 
the orbit. For run A, we terminate the calculation after 12 Myr. 
From the orbital parameters at this point, we estimate that the 
inspiral time due to gravitational wave radiation is of the order 
of a few Myr. Run B terminates as relativistic effects lead to 
SMBH-IMBH coalescence. 



orbit (up to e = 0.97), and GW emission becomes impor- 
tant, leading to further decline of the semimajor axis. After 
3 Myr, GW emission has taken over to rapidly decrease the 
semimajor axis and circularize the orbit to almost e = 0. 
Eventually, the pericentre distance falls below three SMBH 
Schwarzschild radii, and the IMBH is considered to be swal- 
lowed. 



4.3 IMBH inspiral in a flattened cusp 

Recent observations suggest that at least the density pro- 
file of the brightest stars in the Galactic Centre follows 
a more flatt ened cusp than the a = 1.75 profile as sumed 
above (e.g.. iGenzel et all 120031 : ISchodel etail I2007T I. This 
may be the result of mass segregation and the preferen- 
tial depletion of main-se quence stars due to more massive 
stellar mass black holes iBaumgardt. Makino fc Ebisuzakil 
l2004bl : iFreitae. Amaro-Seoane fc Kalogerall2006l'). However^ 
as can b e seen in lBaumgardt. Gualandris fc Portegies Zwartl 
(2006a), this can also be the result of a recent IMBH inspi- 
ral (and thus apply to the whole stellar population in that 
regi on). 

IBaumgardt, Gualandris fc Portegies Zwartl l|2006al ) es- 
timated that replenishment of the c usp after an IMBH in- 
spiral will take at least 100 Myr, while IPortegies Zwart et al.l 
(2006) calculated an IMBH inspiral rate of 1 per 10 Myr into 
the Galactic Centre. This would explain the observed flat- 
tening of the stellar cusp in the centre of the Milky Way, as 
a follow-up IMBH spirals in and scatters away stars before 
the cusp has been replenished. 

In another subsequent model, we study the changes to 
an IMBH inspiral in such a flattened cusp. A second follow- 
up IMBH (run C, A/imbh = 10 3 Mq) finds the density pro- 
file in the inner 0.01 pc to be flattened to oc r -1 ' 2 by the 
preceding two IMBHs, causing the inspiral process to al- 
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Figure 3. Evolution of the inspiralling IMBHs' semimajor axes 
iimbh- The solid, dashed, and dotted lines represent subsequent 
inspiralling IMBHs of mass 10 3 Mq (run A), 3 X 10 3 Mq (run 
B), and 10 3 Mq (run C), respectively. Due to the progressive de- 
pletion of stars in the inner cusp, the inspiral of the subsequent 
IMBHs terminate at larger radii (compare curves for run A and 
C with otherwise identical parameters). The dash-dotted line de- 
picts the 3 X 1O 3 M IMBH in an a = 1.4 cusp (run D). The 
actual inspirals deviate from the theoretical predictions plotted 
in grey once the inner cusp becomes depleted of stars. 



most stop at a distance of 10 3 pc, where the effects of gen- 
eral relativity are still not sufficient to drive a merger with 
the SMBH within 1 Gyr. 

To reproduce the situation of a density profi l e as fo und 
for the centre of the Milky Way bv lGenzel et al.l (|2003l ). we 
also analysed the inspiral of a 3 x 10 3 Mq IMBH in an r~ 1A 
power- law cusp (run D). To follow the dynamical processes 
more accurately, this time we used stellar masses between 
2.5 and 8.5 Mq. According to equation l)17p . for a power law 
with index a < 1.5, the Newtonian decrease of semimajor 
axis will never terminate. 

Fig. [3] shows the evolution of the IMBHs' semimajor 
axes in all four runs. As we have also seen for the previous 
runs, our actual inspiral curve is in good agreement with 
the theoretical prediction for the flattened cusps. The in- 
spiral process is slowed as the IMBH clears out the central 
region. Due to the lower number of stars in the centre of this 
shallower cusp, the decrease of semimajor axis terminates at 
a larger central distance. 



4.4 Detection of IMBH inspirals 

lHansen fc Milosavlievic] (|2003l ) have shown that an IMBH 
of 3 x 10 3 Mq on a circular orbit around Sgr A* at a dis- 
tance between 4 x 10~ 3 and 8 x 10~ 2 pc may be observed 
by detecting the proper motion of Sgr A* if the astromet- 
ric resolution reaches 0.1 mas or better. Fig. U shows the 
amount of time the IMBHs in our simulations spend within 
certain detectability limits. One can see that the IMBH in 
the flattened cusp (run D) spends several Myr in this regime, 
yielding good chance of direct observation in the near future 
if the inspiral rate of IMBHs is large enough. 

At an orbital frequency of 10 -4 Hz, the black hole bi- 
nary enters the LISA detection band. Using the values 
from run B, this corresponds to a semimajor axis of a = 



Figure 4. Detection limits of an astrometric wobble of the radio 
image of Sgr A* as a function of IMBH mass, assuming a circular 
orbit, a SMBH mass of 3 X 10 6 Mq , and a distance to the Galact ic 
Centre of 8.5 kpc (adapted from Hansen & Milosavlievic (2003)). 
The shaded regions show the detection limits in semimajor axis 
a and IMBH mass Mimbh for astrometric resolutions of 10 /ias, 
0.1 mas, and 1 mas, respectively, assuming a monitoring period 
of 10 yr. The dashed and dotted lines indicate coalescence due 
to gravitational wave radiation in 10 7 and 10 s yr. respectively, 
assuming a circular orbit. The motion of the IMBHs from runs A 
to D is represented by solid lines, symbols depict steps of 1 Myr. 
One can see that the 3 x 10 3 Mq IMBH in run D spends 5 Myr 
in a regime where it would be detectable with an astrometric 
resolution of 0.1 mas, and probably a few 10 Myr in the 10 fias 
regime. 



3.25 x 10 -6 pc, and about 2.5 yr until coalescence. Assum- 
ing a distance of R = 1 Gpc from the source, the ampli- 
tude of the GW signal for th e circularized binary is given by 
|Douglas fc Braginskv|[l979l . equation (3.5) and (3.12)) 



h 



128 G 2 MsMBH-AflMBH 

c 4 a7? 



3 x 10" 



(18) 



which is a factor of 3 abov e the noise limit 
(|Seto. Kawamura fc Nakamural 1200 ll ) . Assuming a density 
of SMBHs with masses 1O 6 M < M S mbh < 1O 7 M of 
order dAVdlog M S mbh = 3 x 10 6 Gpc" 3 (lAller fc Richstonei 
2002; iGraham fc Driver! 120071 ) and thus a total number of 
A^smbh ~ 1-3 x 10 7 within a distance of 1 Gpc, as well as 
an IMBH inspiral rate of 1 per 25 Myr, we get a 73 per 
cent chance to find a signal at a given point in time. For 
an observation period of three years, this value rises to 94 
per cent, yielding a good chance to observe IMBH inspirals 
with LISA. Furthermore, SMBHs in other mass ranges will 
contribute to the event rate. 



4.5 HVS ejection 

During the inspiral process, a number of HVSs stars are 
ejected from the system by close encounters with th e IMBH. 
iBaumgardt. Gualandris fc Portegies Zwartl l|2006al ) derived 
relatively short HVS burst intervals of only a few Myr for 
IMBH inspirals in an a = 1.75 power-law cusp. However, 
if the IMBH inspiral terminates at a larger distance (as we 
have seen above for a shallower cusp as found in our Galaxy), 
the ejection rate of HVSs will be lower, and HVSs will be 
created over a longer time interval. For example, the HVS 
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Figure 5. Distribution of galactocontric velocities and distances 
of HVSs created in our simulations. We consider all ejected 
stars in the distance range 10 kpc < r < 120 kpc with velocities 
v > +275 kms" 1 . Line-styles as in Fig. [3] lrr error bars of the 
theoretical distribution are shown in grey and are determined by 
bootstrapping. The filled bars show the d istribution of all B-ty pe 
HVS from the Sloan Digital Sky Survey in lBrown et al.l J2007al fbh. 
Theory would expect 2-3 HVSs with velocities > 1000 kms -1 
which have not been found yet. The low-velocity bin may include 
stars ejected from open clusters. 



ejection rate in run C reaches 4-7 Myr -1 during the first 
5 Myr of HVS ejection, but then retains an average rate of 
3Myr _1 throughout the rest of the simulation. This could 
explain the abse nce of bursts in the observed distribution of 
Galactic HVSs (|Brown et al.ll2007bl ). 

To compare the results of our simulations with the 
HVSs found so far, we followed the trajectories of all ejected 
stars with velocities v > 275 kms -1 at distances between 
10 and 120 kpc from the G alactic Centre. For this range, 
iBrown et al.l (|2006bl . l2007al lbh made a complete survey of 
B-type stars across 12 per cent of the sky. We used a 
Runge-Ku t ta int egrator using the Galactic potential from 
IPaczvnskil l|l99dh adding the broken power-la w density pro- 
file d etermined for the Galactic Centre by iGenzel et al.l 
[|2003r i in the inner 40pc, w here the latter exceeds the den- 
sity of the lPaczvnskH (ll990T ) model. 

We find no significant dependence of the HVS velocities 
produced on the IMBH mass and the density profile of the 
cusp. The distribution of HVSs as they should be observable 
in the Galactic halo is shown in Fig. [5] Our results show that 
there should be a population of HVSs with very high veloc- 
ities (v > 1000 kms" 1 ) which has not been observed so far. 
The observed excess in low-velocity HVSs (v < 500 kms -1 ) 
may be explained by stars ejected from young massive star 
clu sters as a result of binary-binary or IMBH col lisions 
(cf. iGvaramadze. Gualandris fc Portegies Zwarj (|2007f ) and 
references therein). Hence, from the small amount of 
data currently available, the IMBH HVS ejection model 
can neither be excluded nor concluded for our Galaxy. 
These results are consistent wit h the simulations done by 
ISesana. Haardt fc Madaul (120071 ). 

The distribution of central di stances of the HVSs s eems 
compatible with the observations. IBrown et all (|2007af l find 
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Figure 6. IMBH encounter distance for ejected HVSs in units of 
the tidal radius, plotted against the HVS velocity at a distance 
of 10 kpc from the Galactic Centre. The filled squares, crosses, 
asterisks, and open circles represent the data from run A, B, 
C, and D, respectively. Most HVSs created passed the IMBH 
at a distance of more than several tidal radii and are thus not 
affected by tidal effects. Only among the very fast moving HVSs 
(v > 1000 kms -1 ), a significant proportion has passed the IMBH 
within a few tidal radii. Stars passing the IMBH below iJ t id will 
get tidally disrupted. 



a sligh t ani sotropy in the spatial distribution of HVSs; lLevinl 
l|200rj| and lSesana. Haardt fc Madaul (|200rj ') argue that the 
IMBH's orbital plane gives preferred directions for HVS ejec- 
tion. However, our simulations show that the orbital plane 
changes rapidly over the inspiral time of a few Myr. Thus, 
the distribution of HVS directions may be fairly random, 
even if caused by only a few IMBH inspirals. 

4.6 HVS rotation 

HVSs are created in the Galactic Centre as a result of close 
encounters with massive black holes. If these encounters 
are close enough, tidal effects become important. One ef- 
fect would be that, as a result of the encounter, rotation is 
induced in the star. 

In the case of stellar binaries disrupted by the central 
SMBH (|Hilldll988h . a binary star with semimajor axis a and 
masses Mi , M% will get disrupted at 



iitid = a ^/A./smbh/ (Mi + M 2 ) (19) 
l|Yu fc Trem ain c 2003). In terms of tidal radii of the member 



stars, iitid = \/A/smbh/M,, R t , this is R^/Rtid = 0.8 a/ R* 
for an equal-mass binary. Thus, if a binary's semimajor axis 
is at least a few times the radius of each member, a HVS 
created by binary disruption passes the SMBH no closer 
than a few tidal radii, meaning that it would not feel a strong 
tidal force due to the SMBH. 

In the case of HVSs created by IMBH encounters, the 
tidal disruption radius of a B-type star with radius i?„ = 
1.5 R® and mass M, = 3 M Q is R tid = 10.4 R Q (15i? Q ), 
assuming a IMBH mass of Mimbh = 10 3 M Q (3 x 10 3 M ). 

Fig.|6]shows the IMBH encounter distance of the HVSs 
created in our simulations. Only among the very fast moving 
HVSs (v > 1000 kms -1 ), a significant proportion has passed 
the IMBH within a few tidal radii. 
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Figure 7. Induced rotation of a star due to a close passage of 
an IMBH. The curves are labelled by the minimum encounter 
distance (t = 0) in terms of the star's tidal radius R tl d- Significant 
rotation is only induced if the minimum distance is r m j n < ZRtid ■ 

In the following, we calculate the amount of rotation 
induced by a massive black hole passage, depending on the 
minimum distance. Detection of rapid rotation of any ob- 
served HVS could be evidence for the presence of IMBHs in 
the Galactic Centre. 

In order to test up to which encounter distance rotation 
is induced in a star, we m ade calculations with a mod ified 
version of the SPH code of iNakasato fc Nomotol (|2003h . We 
simulated hyperbolic encounters of a 3 Mq star with a radius 
of R* = 1.5 Rq and a 3 x 10 3 Af© IMBH. Fig. shows that 
encounters with a minimum distance of less than two tidal 
radii from the IMBH cause the passing star to rotate signif- 
icantly. An encounter within 1.3 tidal radii induces rotation 
with a period of two hours. 

A significant rotation should be visible for those stars 
passing the IMBH within two tidal radii, which in our calcu- 
lations are only 3 per cent of the total number of HVSs cre- 
ated. Hence, only among very fast HVSs (v > 1000 kms -1 ) 
we can expect to find a significant proportion of rotating 
stars as evidence for an IMBH in the Galactic Centre. 

4.7 Ejection of hypervelocity binary stars 

iLu. Yu fc Linl (|2007T ) showed that binary stars can be ejected 
from the Galactic Centre as a whole through interactions 
with an IMBH. Using equation (|19p with Mimbh = 3 x 
10 3 M Q and Mi = M 2 = 3M , we get « 8a as the 

minimum distance from the IMBH which a binary star with 
semimajor axis a could survive before being tidally dis- 
rupted. As one can see in Fig. [6] all HVSs created in our 
models have encounter distances below 400 Rtid, or 30 au. 
Hence, to survive these encounters, the binary stars must 
have a separation below 3.75 au. 

To survive the tidal forces at a distance of 0.05 pc from 
a 3 x 10 6 Mq SMBH, a 6 Mq binary must have a semimajor 
axis below 130 au, giving an upper limit for the separation 
of stellar binaries in the Galactic Centre. So far, not much 
is known about the binary fraction /bin and the separation 
distribution in the Galactic Centre. As a rough estimate, we 
can use the period distribution of G-dwarf binaries in the 



solar neighbourhood l|Duauennov fc Mavor|[l99lh . We find 
that 37 per cent of all binaries with semimajor axis below 
130 au have a separation below 3.75 au and may eventually 
be ejected as hypervelocity binary stars (HVBSs). 

Binaries ejected tend to be among the low-velocity 
HVSs. If we consider the actual distribution of encounter 
distances as shown in Fig. [5] and also the different 'life- 
times' of the HVSs, i.e. the amount of time a star spends in 
the HVS regime with velocity v > 275 kms -1 and distance 
to the Galactic Centre lOkpc < r < 120 kpc, we find that 
the actual fraction of HVBSs is only 0.13 x /bin, i.e. no more 
than a few per cent. 



5 DISCUSSION 

We have developed a novel post-Newtonian, high-precision 
integrator bhint for stellar systems containing a SMBH. The 
algorithm makes use of the fact that the Keplerian orbits 
in such a potential can be calculated directly, and are only 
weakly perturbed. We used this method to follow the inspiral 
of IMBHs into the Galactic Centre as a result of dynamical 
friction. We find that on their way to the centre, IMBHs 
deplete the central cusp of stars, leaving behind a shallower 
density profile compatible with the one recently observed in 
the centre of the Milky Way. Due to the decreased stellar 
density, IMBHs may get stuck at a distance of a milliparsec 
from the SMBH for a few 10 Myr. If the recently detected 
HVSs have been ejected by one or several IMBHs, it is likely 
that an IMBH still resides close to Sgr A*. 

Such an IMBH may be observable by detecting the 
proper motion of the radio source Sgr A* with upcom- 
ing high-resolution telescopes. Our results suggest that an 
IMBH in the Galactic Centre may spend several Myr in the 
regime where it is detectable with an astrometric resolution 
of 0.1 mas, yielding good chance of direct observation in the 
near future. 

Once an IMBH comes close enough to the centre, gravi- 
tational wave emission will take over and lead to the merger 
with the central SMBH. These mergers are among the most 
powerful sources of gr avitational waves for the planned LISA 
mission (|Thorneil987l ) and may well be detected if the IMBH 
inspiral rate is of the order of 1 per 25 Myr. 

Our results show that there should be a population of 
HVSs with very high velocities (v > 1000 kms -1 ) which 
h as not been observed so far, con sistent with the findings 
of lSesana. Haardt fc Madaul(|2007T l. However, as the sample 
size of HVSs observed is relatively small, IMBHs cannot be 
excluded as ejection mechanism. Moreover, HVSs ejected by 
IMBHs are hard to distinguish from those ejected by other 
mechanisms such as binary disruptions as we do not expect 
significant rotation in either case, but possible observations 
of significant rotation of HVSs in the high velocity tail may 
favour the IMBH model. 



ACKNOWLEDGEMENTS 

We are grateful to Douglas C. Heggie and Ulrich Heber for 
useful discussions, and to Ingo Berentzen for his help with 
the implementation of the post-Newtonian extension. We 
also thank Naohito Nakasato for the help with his SPH code. 



Tracing IMBHs in the Galactic Centre 9 



This work was supported by the DFG Priority Program 1177 
'Witnesses of Cosmic History: Formation and Evolution of 
Black Holes, Galaxies and Their Environment'. 



REFERENCES 

Aarseth S. J., 1985, in Brackbill J. U., Cohen B. I., eds, 
Multiple Time Scales. Academic Press, New York, p. 377 
Aarseth S. J., 2007, MNRAS, 378, 285 
Aller M. C, Richstone D., 2002, A J, 124, 3035 
Bahcall J. N., Wolf R. A., 1976, ApJ, 209, 214 
Baumgardt H., Makino J., Ebisuzaki T., 2004a, ApJ, 613, 
1133 

Baumgardt H., Makino J., Ebisuzaki T., 2004b, ApJ, 613, 
1143 

Baumgardt H., Gualandris A., Portegies Zwart S., 2006a, 
MNRAS, 372, 174 

Baumgardt H., Hopman C, Portegies Zwart S., Makino J., 
2006b, MNRAS, 372, 467 

Bender P. et al., 1998, LISA. Laser Interferometer Space 
Antenna for the detection and observation of gravita- 
tional waves. Pre-Phase A report, Tech. rep., Max-Planck- 
Institut fur Quantenoptik, Garching 

Binney J., Tremaine S., 1987, Galactic dynamics. Princeton 
Univ. Press, Princeton, NJ 



Blanchet L., 2006, Living Reviews in Relativity, 9, 4 



Brown W. R., Geller M. J., 


Kenyon 


S. 


J.. 


Kurtz 


M. 


J.. 


2005, ApJ, 622, L33 














Brown W. R., Geller M. J., 


Kenyon 


S. 


J., 


Kurtz 


M. 


J., 


2006a, ApJ, 640, L35 














Brown W. R., Geller M. J., 


Kenyon 


S. 


J., 


Kurtz 


M. 


J., 


2006b, ApJ, 647, 303 














Brown W. R., Geller M. J., 


Kenyon 


s. 


J., 


Kurtz 


M. 


J., 


Bromley B. C, 2007a, ApJ, 660, 311 












Brown W. R., Geller M. J., 


Kenyon 


s. 


J., 


Kurtz 


M. 


J., 


Bromley B.C., 2007b, ApJ, 


in press (arXiv:0709.1471 


) 



Chatterjee P., Hernquist L., Loeb A., 2002, ApJ, 572, 371 
Douglas D. H., Braginsky V. G., 1979, in Hawking S. W., 
Israel W., eds, General Relativity: An Einstein Centenary 
Survey. Cambridge Univ. Press, Cambridge, p. 90 



Duquennoy A., Mayor M., 1991, A&A, 248, 485 
Edelmann H., Napiwotzki R., Heber U., Christlieb N., 

Reimers D., 2005, ApJ, 634, L181 
Freitag M., Amaro-Seoane P., Kalogera V., 2006, ApJ, 649, 

91 

Fujii M., Iwasawa M., Funato Y., Makino J., 2007, PASJ, 

accepted (astro-ph/0706.2059) 
Gebhardt K. et al., 2000, ApJ, 539, L13 
Genzel R. et al., 2003, ApJ, 594, 812 
Ghez A. M., Salim S., Hornstein S. D., Tanner A., Lu J. R., 

Morris M., Becklin E. E., Duchene G., 2005, ApJ, 620, 744 
Graham A. W., Driver S. P., 2007, MNRAS, 380, L15 
Giirkan M. A., Freitag M., Rasio F. A., 2004, ApJ, 604, 

632 

Gvaramadze V. V., Gualandris A., Portegies Zwart S., 

2007, preprint ( astro-ph/0702735[ ) 
Hansen B. M. S., Milosavljevic M., 2003, ApJ, 593, L77 
Hills J. G., 1988, Nat, 331, 687 

Hirsch H. A., Heber U., O'Toole S. J., Bresolin F., 2005, 

A&A, 444, L61 
Hopman C, Alexander T., 2006, ApJ, 645, 1152 



Hopman C, Portegies Zwart S. F., Alexander T., 2004, 

ApJ, 604, L101 
Kokubo E., Yoshinaga K., Makino J., 1998, MNRAS, 297, 

1067 

Kupi G., Amaro-Seoane P., Spurzem R., 2006, MNRAS, 

371, L45 
Levin Y., 2006, ApJ, 653, 1203 

Levin Y., Wu A., Thommes E., 2005, ApJ, 635, 341 

Lu Y., Yu Q., Lin D. N. C, 2007, ApJ, 666, L89 

Magorrian J. et al., 1998, AJ, 115, 2285 

Makino J., Aarseth S. J., 1992, PASJ, 44, 141 

Makino J., Fukushige T., Koga M., Namura K., 2003, 

PASJ, 55, 1163 
Matsubayashi T., Makino J., Ebisuzaki T., 2007, ApJ, 656, 

879 

Mikkola S., Tanikawa K., 1999, Celestial Mechanics and 

Dynamical Astronomy, 74, 287 
Miller M. C, Colbert E. J. M., 2004, International Journal 

of Modern Physics D, 13, 1 
Miralda-Escude J., Gould A., 2000, ApJ, 545, 847 
Nakasato N., Nomoto K., 2003, ApJ, 588, 842 
O'Leary R. M., Loeb A., 2008, MNRAS, 383, 86 
Paczynski B., 1990, ApJ, 348, 485 

Perets H. B., Hopman C, Alexander T., 2007, ApJ, 656, 
709 

Peters P. C, 1964, Physical Review, 136, 1224 

Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., 
McMillan S. L. W., 2004, Nat, 428, 724 

Portegies Zwart S. F., Baumgardt H., McMillan S. L. W., 
Makino J., Hut P., Ebisuzaki T., 2006, ApJ, 641, 319 

Preto M., Merritt D., Spurzem R., 2004, ApJ, 613, L109 

Saha P., Tremaine S., 1994, AJ, 108, 1962 

Schodel R. et al., 2007, A&A, 469, 125 

Schodel R., Ott T., Genzel R., Eckart A., Mouawad N., 
Alexander T., 2003, ApJ, 596, 1015 

Sesana A., Haardt F., Madau P., 2006, ApJ, 651, 392 

Sesana A., Haardt F., Madau P., 2007, MNRAS, L53+ 

Seto N., Kawamura S., Nakamura T., 2001, Phys. Rev. 
Lett., 87, 221103 

Thorne K. S., 1987, Gravitational radiation (Three hun- 
dred years of gravitation), 330-458 

Wisdom J., Holman M., 1991, AJ, 102, 1528 

Yu Q., Tremaine S., 2003, ApJ, 599, 1129 



